home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Turnbull China Bikeride
/
Turnbull China Bikeride - Disc 1.iso
/
ARGONET
/
PD
/
MATHS
/
RLAB
/
RLAB125.ZIP
/
!RLaB
/
examples
/
house
< prev
next >
Wrap
Text File
|
1994-10-22
|
2KB
|
116 lines
//
// From MATRIX Computations, G.H. Golub, C.F. Van Loan
//
//
// house_v(): Given an N-vector V, generate an n-vector V
// with V[1] = 1, such that (eye(n,n) - 2*(V*V')/(V'*V))*X
// is zero in all but the 1st component.
//
static (sign)
house_v = function(v)
{
local(v)
n = max( size(v) );
u = norm(v, "2");
if(u != 0)
{
b = v[1] + sign(v[1])*u;
if(n > 1)
{
v[2:n] = v[2:n]/b;
}
}
v[1] = 1;
return v;
};
//
// house_row(): Given an MxN matrix A and a non-zero M-vector V
// with V[1] = 1, the following algorithm overwrites A with
// P*A, where P = eye(m,m) - 2*(V*V')/(V'*V)
//
house_row = function(A, v)
{
local(A)
b = -2/(v'*v);
w = b*A'*v;
A = A + v*w';
return A;
};
//
// house_col(): Given an MxN matrix A, and an N-vector V,
// with V[1] = 1, the following algorithm overwrites A with A*P
// where P = eye(N,N) - 2*(V*V')/(V'*V)
//
house_col = function(A, v)
{
local(A)
b = -2/(v'*v);
w = b*A*v;
a = A + w*v';
return A;
};
//
// Given A, with M >= N, the following function finds Householder
// matrices H1,...Hn, such that if Q = H1*...Hn, then Q'*A = R is
// upper triangular.
// House_qr returns a MxN matrix, with the upper triangular part
// containing [R]
house_qr = function ( A )
{
local (A)
m = A.nr; n = A.nc;
v = zeros(m,1);
q = eye (m, m);
for(j in 1:n)
{
// Generate the Householder vector
v[j:m] = house_v( A[j:m;j] );
// Apply the Householder reflector to A
A[j:m;j:n] = house_row( A[j:m;j:n], v[j:m] );
// Create Q
if(j < m)
{
q = P ([ zeros (j-1,1); 1; v[j+1:m] ]) * q;
}
}
return << q = q'; r = A >>;
};
//
// Generate P matrix
//
P = function ( V )
{
m = max( size(V) );
return [ eye(m,m) - 2*(V*V')./(V'*V) ];
};
sign = function ( X )
{
if (X >= 0)
{
return 1;
else
return -1;
}
};